Signal Decay

In this chapter, we use simulated data to show how BOLD signal decays and how we can glean useful information from that decay rate.

import os

import imageio
import matplotlib.pyplot as plt
import nibabel as nib
import numpy as np
import seaborn as sns
from IPython import display
from matplotlib.animation import FuncAnimation
from myst_nb import glue
from nilearn import image, plotting
from nilearn.glm import first_level
from repo2data.repo2data import Repo2Data
from scipy import signal

from book_utils import predict_bold_signal

sns.set_style("whitegrid")
plt.rcParams.update({
    "text.usetex": True,
    "font.family": "sans-serif",
    "font.sans-serif": ["Helvetica"]
})

# Install the data if running locally, or point to cached data if running on neurolibre
DATA_REQ_FILE = os.path.join("../binder/data_requirement.json")

# Download data
repo2data = Repo2Data(DATA_REQ_FILE)
data_path = repo2data.install()
data_path = os.path.abspath(os.path.join(data_path[0], "data"))

out_dir = os.path.join(data_path, "signal-decay")
os.makedirs(out_dir, exist_ok=True)
---- repo2data starting ----
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/repo2data
Config from file :
../binder/data_requirement.json
Destination:
./../data/multi-echo-data-analysis

Info : ./../data/multi-echo-data-analysis already downloaded
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/nilearn/datasets/__init__.py:96: FutureWarning: Fetchers from the nilearn.datasets module will be updated in version 0.9 to return python strings instead of bytes and Pandas dataframes instead of Numpy arrays.
  "Numpy arrays.", FutureWarning)

Signal decays as echo time increases

func_dir = os.path.join(data_path, "sub-04570/func/")
data_files = [
    os.path.join(func_dir, "sub-04570_task-rest_echo-1_space-scanner_desc-partialPreproc_bold.nii.gz"),
    os.path.join(func_dir, "sub-04570_task-rest_echo-2_space-scanner_desc-partialPreproc_bold.nii.gz"),
    os.path.join(func_dir, "sub-04570_task-rest_echo-3_space-scanner_desc-partialPreproc_bold.nii.gz"),
    os.path.join(func_dir, "sub-04570_task-rest_echo-4_space-scanner_desc-partialPreproc_bold.nii.gz"),
]
echo_times = [12., 28., 44., 60.]

img = image.index_img(data_files[0], 0)
data = img.get_fdata()
vmax = np.max(data)
idx = np.vstack(np.where(data > 1000))

min_x, min_y, min_z = np.min(idx, axis=1)
max_x, max_y, max_z = np.max(idx, axis=1)

imgs = []
for f in data_files:
    img = image.index_img(f, 0)
    data = img.get_fdata()
    data = data[min_x:max_x, min_y:max_y, min_z:max_z]
    img = nib.Nifti1Image(data, img.affine, img.header)
    imgs.append(img)
plt.style.use("dark_background")

fig, axes = plt.subplots(figsize=(26, 4), ncols=len(data_files))
for i_echo, img in enumerate(imgs):
    te = echo_times[i_echo]
    if i_echo == 0:
        data = img.get_fdata()
        vmax = np.max(data)

    plotting.plot_epi(
        img,
        cut_coords=[0],
        display_mode="z",
        annotate=False,
        draw_cross=False,
        axes=axes[i_echo],
        figure=fig,
        vmax=vmax,
        cmap="gray",
    )
    axes[i_echo].set_title(f"TE={te}ms", fontsize=20, pad=0)

glue("fig_brain_decay", fig, display=False)

# Reset the style
plt.style.use("default")
---------------------------------------------------------------------------
FileNotFoundError                         Traceback (most recent call last)
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/texmanager.py in _run_checked_subprocess(self, command, tex, cwd)
    253                 command, cwd=cwd if cwd is not None else self.texcache,
--> 254                 stderr=subprocess.STDOUT)
    255         except FileNotFoundError as exc:

/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/subprocess.py in check_output(timeout, *popenargs, **kwargs)
    410     return run(*popenargs, stdout=PIPE, timeout=timeout, check=True,
--> 411                **kwargs).stdout
    412 

/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/subprocess.py in run(input, capture_output, timeout, check, *popenargs, **kwargs)
    487 
--> 488     with Popen(*popenargs, **kwargs) as process:
    489         try:

/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/subprocess.py in __init__(self, args, bufsize, executable, stdin, stdout, stderr, preexec_fn, close_fds, shell, cwd, env, universal_newlines, startupinfo, creationflags, restore_signals, start_new_session, pass_fds, encoding, errors, text)
    799                                 errread, errwrite,
--> 800                                 restore_signals, start_new_session)
    801         except:

/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/subprocess.py in _execute_child(self, args, executable, preexec_fn, close_fds, pass_fds, cwd, env, startupinfo, creationflags, shell, p2cread, p2cwrite, c2pread, c2pwrite, errread, errwrite, restore_signals, start_new_session)
   1550                             err_msg += ': ' + repr(err_filename)
-> 1551                     raise child_exception_type(errno_num, err_msg, err_filename)
   1552                 raise child_exception_type(err_msg)

FileNotFoundError: [Errno 2] No such file or directory: 'latex': 'latex'

The above exception was the direct cause of the following exception:

RuntimeError                              Traceback (most recent call last)
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/IPython/core/formatters.py in __call__(self, obj)
    339                 pass
    340             else:
--> 341                 return printer(obj)
    342             # Finally look for special method names
    343             method = get_real_method(obj, self.print_method)

/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/IPython/core/pylabtools.py in print_figure(fig, fmt, bbox_inches, base64, **kwargs)
    149         FigureCanvasBase(fig)
    150 
--> 151     fig.canvas.print_figure(bytes_io, **kw)
    152     data = bytes_io.getvalue()
    153     if fmt == 'svg':

/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/backend_bases.py in print_figure(self, filename, dpi, facecolor, edgecolor, orientation, format, bbox_inches, pad_inches, bbox_extra_artists, backend, **kwargs)
   2228                        else suppress())
   2229                 with ctx:
-> 2230                     self.figure.draw(renderer)
   2231 
   2232             if bbox_inches:

/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/artist.py in draw_wrapper(artist, renderer, *args, **kwargs)
     72     @wraps(draw)
     73     def draw_wrapper(artist, renderer, *args, **kwargs):
---> 74         result = draw(artist, renderer, *args, **kwargs)
     75         if renderer._rasterizing:
     76             renderer.stop_rasterizing()

/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/artist.py in draw_wrapper(artist, renderer, *args, **kwargs)
     49                 renderer.start_filter()
     50 
---> 51             return draw(artist, renderer, *args, **kwargs)
     52         finally:
     53             if artist.get_agg_filter() is not None:

/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/figure.py in draw(self, renderer)
   2779             self.patch.draw(renderer)
   2780             mimage._draw_list_compositing_images(
-> 2781                 renderer, self, artists, self.suppressComposite)
   2782 
   2783             for sfig in self.subfigs:

/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/image.py in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
    130     if not_composite or not has_images:
    131         for a in artists:
--> 132             a.draw(renderer)
    133     else:
    134         # Composite any adjacent images together

/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/artist.py in draw_wrapper(artist, renderer, *args, **kwargs)
     49                 renderer.start_filter()
     50 
---> 51             return draw(artist, renderer, *args, **kwargs)
     52         finally:
     53             if artist.get_agg_filter() is not None:

/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/_api/deprecation.py in wrapper(*inner_args, **inner_kwargs)
    429                          else deprecation_addendum,
    430                 **kwargs)
--> 431         return func(*inner_args, **inner_kwargs)
    432 
    433     return wrapper

/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/axes/_base.py in draw(self, renderer, inframe)
   2879                 artists.remove(spine)
   2880 
-> 2881         self._update_title_position(renderer)
   2882 
   2883         if not self.axison or inframe:

/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/axes/_base.py in _update_title_position(self, renderer)
   2810                 if (ax.xaxis.get_ticks_position() in ['top', 'unknown']
   2811                         or ax.xaxis.get_label_position() == 'top'):
-> 2812                     bb = ax.xaxis.get_tightbbox(renderer)
   2813                 else:
   2814                     bb = ax.get_window_extent(renderer)

/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/axis.py in get_tightbbox(self, renderer, for_layout_only)
   1081         ticks_to_draw = self._update_ticks()
   1082 
-> 1083         self._update_label_position(renderer)
   1084 
   1085         # go back to just this axis's tick labels

/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/axis.py in _update_label_position(self, renderer)
   2078         # get bounding boxes for this axis and any siblings
   2079         # that have been set by `fig.align_xlabels()`
-> 2080         bboxes, bboxes2 = self._get_tick_boxes_siblings(renderer=renderer)
   2081 
   2082         x, y = self.label.get_position()

/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/axis.py in _get_tick_boxes_siblings(self, renderer)
   1866             axis = ax._get_axis_map()[axis_name]
   1867             ticks_to_draw = axis._update_ticks()
-> 1868             tlb, tlb2 = axis._get_tick_bboxes(ticks_to_draw, renderer)
   1869             bboxes.extend(tlb)
   1870             bboxes2.extend(tlb2)

/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/axis.py in _get_tick_bboxes(self, ticks, renderer)
   1062         """Return lists of bboxes for ticks' label1's and label2's."""
   1063         return ([tick.label1.get_window_extent(renderer)
-> 1064                  for tick in ticks if tick.label1.get_visible()],
   1065                 [tick.label2.get_window_extent(renderer)
   1066                  for tick in ticks if tick.label2.get_visible()])

/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/axis.py in <listcomp>(.0)
   1062         """Return lists of bboxes for ticks' label1's and label2's."""
   1063         return ([tick.label1.get_window_extent(renderer)
-> 1064                  for tick in ticks if tick.label1.get_visible()],
   1065                 [tick.label2.get_window_extent(renderer)
   1066                  for tick in ticks if tick.label2.get_visible()])

/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/text.py in get_window_extent(self, renderer, dpi)
    901 
    902         with cbook._setattr_cm(self.figure, dpi=dpi):
--> 903             bbox, info, descent = self._get_layout(self._renderer)
    904             x, y = self.get_unitless_position()
    905             x, y = self.get_transform().transform((x, y))

/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/text.py in _get_layout(self, renderer)
    306         _, lp_h, lp_d = renderer.get_text_width_height_descent(
    307             "lp", self._fontproperties,
--> 308             ismath="TeX" if self.get_usetex() else False)
    309         min_dy = (lp_h - lp_d) * self._linespacing
    310 

/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/backends/backend_agg.py in get_text_width_height_descent(self, s, prop, ismath)
    228             fontsize = prop.get_size_in_points()
    229             w, h, d = texmanager.get_text_width_height_descent(
--> 230                 s, fontsize, renderer=self)
    231             return w, h, d
    232 

/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/texmanager.py in get_text_width_height_descent(self, tex, fontsize, renderer)
    397         else:
    398             # use dviread.
--> 399             dvifile = self.make_dvi(tex, fontsize)
    400             with dviread.Dvi(dvifile, 72 * dpi_fraction) as dvi:
    401                 page, = dvi

/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/texmanager.py in make_dvi(self, tex, fontsize)
    291                 self._run_checked_subprocess(
    292                     ["latex", "-interaction=nonstopmode", "--halt-on-error",
--> 293                      texfile], tex, cwd=tmpdir)
    294                 (Path(tmpdir) / Path(dvifile).name).replace(dvifile)
    295         return dvifile

/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/texmanager.py in _run_checked_subprocess(self, command, tex, cwd)
    256             raise RuntimeError(
    257                 'Failed to process string with tex because {} could not be '
--> 258                 'found'.format(command[0])) from exc
    259         except subprocess.CalledProcessError as exc:
    260             raise RuntimeError(

RuntimeError: Failed to process string with tex because latex could not be found
---------------------------------------------------------------------------
FileNotFoundError                         Traceback (most recent call last)
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/texmanager.py in _run_checked_subprocess(self, command, tex, cwd)
    253                 command, cwd=cwd if cwd is not None else self.texcache,
--> 254                 stderr=subprocess.STDOUT)
    255         except FileNotFoundError as exc:

/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/subprocess.py in check_output(timeout, *popenargs, **kwargs)
    410     return run(*popenargs, stdout=PIPE, timeout=timeout, check=True,
--> 411                **kwargs).stdout
    412 

/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/subprocess.py in run(input, capture_output, timeout, check, *popenargs, **kwargs)
    487 
--> 488     with Popen(*popenargs, **kwargs) as process:
    489         try:

/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/subprocess.py in __init__(self, args, bufsize, executable, stdin, stdout, stderr, preexec_fn, close_fds, shell, cwd, env, universal_newlines, startupinfo, creationflags, restore_signals, start_new_session, pass_fds, encoding, errors, text)
    799                                 errread, errwrite,
--> 800                                 restore_signals, start_new_session)
    801         except:

/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/subprocess.py in _execute_child(self, args, executable, preexec_fn, close_fds, pass_fds, cwd, env, startupinfo, creationflags, shell, p2cread, p2cwrite, c2pread, c2pwrite, errread, errwrite, restore_signals, start_new_session)
   1550                             err_msg += ': ' + repr(err_filename)
-> 1551                     raise child_exception_type(errno_num, err_msg, err_filename)
   1552                 raise child_exception_type(err_msg)

FileNotFoundError: [Errno 2] No such file or directory: 'latex': 'latex'

The above exception was the direct cause of the following exception:

RuntimeError                              Traceback (most recent call last)
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/IPython/core/formatters.py in __call__(self, obj)
    339                 pass
    340             else:
--> 341                 return printer(obj)
    342             # Finally look for special method names
    343             method = get_real_method(obj, self.print_method)

/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/IPython/core/pylabtools.py in print_figure(fig, fmt, bbox_inches, base64, **kwargs)
    149         FigureCanvasBase(fig)
    150 
--> 151     fig.canvas.print_figure(bytes_io, **kw)
    152     data = bytes_io.getvalue()
    153     if fmt == 'svg':

/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/backend_bases.py in print_figure(self, filename, dpi, facecolor, edgecolor, orientation, format, bbox_inches, pad_inches, bbox_extra_artists, backend, **kwargs)
   2228                        else suppress())
   2229                 with ctx:
-> 2230                     self.figure.draw(renderer)
   2231 
   2232             if bbox_inches:

/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/artist.py in draw_wrapper(artist, renderer, *args, **kwargs)
     72     @wraps(draw)
     73     def draw_wrapper(artist, renderer, *args, **kwargs):
---> 74         result = draw(artist, renderer, *args, **kwargs)
     75         if renderer._rasterizing:
     76             renderer.stop_rasterizing()

/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/artist.py in draw_wrapper(artist, renderer, *args, **kwargs)
     49                 renderer.start_filter()
     50 
---> 51             return draw(artist, renderer, *args, **kwargs)
     52         finally:
     53             if artist.get_agg_filter() is not None:

/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/figure.py in draw(self, renderer)
   2779             self.patch.draw(renderer)
   2780             mimage._draw_list_compositing_images(
-> 2781                 renderer, self, artists, self.suppressComposite)
   2782 
   2783             for sfig in self.subfigs:

/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/image.py in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
    130     if not_composite or not has_images:
    131         for a in artists:
--> 132             a.draw(renderer)
    133     else:
    134         # Composite any adjacent images together

/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/artist.py in draw_wrapper(artist, renderer, *args, **kwargs)
     49                 renderer.start_filter()
     50 
---> 51             return draw(artist, renderer, *args, **kwargs)
     52         finally:
     53             if artist.get_agg_filter() is not None:

/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/_api/deprecation.py in wrapper(*inner_args, **inner_kwargs)
    429                          else deprecation_addendum,
    430                 **kwargs)
--> 431         return func(*inner_args, **inner_kwargs)
    432 
    433     return wrapper

/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/axes/_base.py in draw(self, renderer, inframe)
   2879                 artists.remove(spine)
   2880 
-> 2881         self._update_title_position(renderer)
   2882 
   2883         if not self.axison or inframe:

/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/axes/_base.py in _update_title_position(self, renderer)
   2810                 if (ax.xaxis.get_ticks_position() in ['top', 'unknown']
   2811                         or ax.xaxis.get_label_position() == 'top'):
-> 2812                     bb = ax.xaxis.get_tightbbox(renderer)
   2813                 else:
   2814                     bb = ax.get_window_extent(renderer)

/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/axis.py in get_tightbbox(self, renderer, for_layout_only)
   1081         ticks_to_draw = self._update_ticks()
   1082 
-> 1083         self._update_label_position(renderer)
   1084 
   1085         # go back to just this axis's tick labels

/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/axis.py in _update_label_position(self, renderer)
   2078         # get bounding boxes for this axis and any siblings
   2079         # that have been set by `fig.align_xlabels()`
-> 2080         bboxes, bboxes2 = self._get_tick_boxes_siblings(renderer=renderer)
   2081 
   2082         x, y = self.label.get_position()

/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/axis.py in _get_tick_boxes_siblings(self, renderer)
   1866             axis = ax._get_axis_map()[axis_name]
   1867             ticks_to_draw = axis._update_ticks()
-> 1868             tlb, tlb2 = axis._get_tick_bboxes(ticks_to_draw, renderer)
   1869             bboxes.extend(tlb)
   1870             bboxes2.extend(tlb2)

/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/axis.py in _get_tick_bboxes(self, ticks, renderer)
   1062         """Return lists of bboxes for ticks' label1's and label2's."""
   1063         return ([tick.label1.get_window_extent(renderer)
-> 1064                  for tick in ticks if tick.label1.get_visible()],
   1065                 [tick.label2.get_window_extent(renderer)
   1066                  for tick in ticks if tick.label2.get_visible()])

/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/axis.py in <listcomp>(.0)
   1062         """Return lists of bboxes for ticks' label1's and label2's."""
   1063         return ([tick.label1.get_window_extent(renderer)
-> 1064                  for tick in ticks if tick.label1.get_visible()],
   1065                 [tick.label2.get_window_extent(renderer)
   1066                  for tick in ticks if tick.label2.get_visible()])

/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/text.py in get_window_extent(self, renderer, dpi)
    901 
    902         with cbook._setattr_cm(self.figure, dpi=dpi):
--> 903             bbox, info, descent = self._get_layout(self._renderer)
    904             x, y = self.get_unitless_position()
    905             x, y = self.get_transform().transform((x, y))

/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/text.py in _get_layout(self, renderer)
    306         _, lp_h, lp_d = renderer.get_text_width_height_descent(
    307             "lp", self._fontproperties,
--> 308             ismath="TeX" if self.get_usetex() else False)
    309         min_dy = (lp_h - lp_d) * self._linespacing
    310 

/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/backends/backend_agg.py in get_text_width_height_descent(self, s, prop, ismath)
    228             fontsize = prop.get_size_in_points()
    229             w, h, d = texmanager.get_text_width_height_descent(
--> 230                 s, fontsize, renderer=self)
    231             return w, h, d
    232 

/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/texmanager.py in get_text_width_height_descent(self, tex, fontsize, renderer)
    397         else:
    398             # use dviread.
--> 399             dvifile = self.make_dvi(tex, fontsize)
    400             with dviread.Dvi(dvifile, 72 * dpi_fraction) as dvi:
    401                 page, = dvi

/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/texmanager.py in make_dvi(self, tex, fontsize)
    291                 self._run_checked_subprocess(
    292                     ["latex", "-interaction=nonstopmode", "--halt-on-error",
--> 293                      texfile], tex, cwd=tmpdir)
    294                 (Path(tmpdir) / Path(dvifile).name).replace(dvifile)
    295         return dvifile

/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/matplotlib/texmanager.py in _run_checked_subprocess(self, command, tex, cwd)
    256             raise RuntimeError(
    257                 'Failed to process string with tex because {} could not be '
--> 258                 'found'.format(command[0])) from exc
    259         except subprocess.CalledProcessError as exc:
    260             raise RuntimeError(

RuntimeError: Failed to process string with tex because latex could not be found
<Figure size 1872x288 with 8 Axes>
<Figure size 1872x288 with 8 Axes>

Fig. 1 Signal decay in the brain.

The impact of \(T_{2}^{*}\) and \(S_{0}\) fluctuations on BOLD signal

# Simulate data
MULTIECHO_TES = np.array([15, 30, 45, 60, 75, 90])
SINGLEECHO_TE = np.array([30])

# For a nice, smooth curve
FULLCURVE_TES = np.arange(0, 101, 1)

n_echoes = len(FULLCURVE_TES)
pal = sns.color_palette('cubehelix', 8)

Simulate \(T_{2}^{*}\) and \(S_{0}\) fluctuations

# Simulate data
# We'll convolve with HRF just for smoothness
hrf = first_level.spm_hrf(1, oversampling=1)

N_VOLS = 21

SCALING_FRACTION = 0.1  # used to scale standard deviation
MEAN_T2S = 35
MEAN_S0 = 16000

# simulate the T2*/S0 time series
# The original time series will be a random time series from a normal distribution,
# convolved with the HRF
ts = np.random.normal(loc=0, scale=1, size=(N_VOLS+20,))
ts = signal.convolve(ts, hrf)[20:N_VOLS+20]
ts *= SCALING_FRACTION / np.std(ts)
ts -= np.mean(ts)

t2s_ts = (ts * MEAN_T2S) + MEAN_T2S
s0_ts = (ts * MEAN_S0) + MEAN_S0

Plot BOLD signal decay for a standard single-echo scan

fullcurve_signal = predict_bold_signal(FULLCURVE_TES, [MEAN_S0], [30])

fig, ax = plt.subplots(figsize=(14, 7.5))
ax.plot(
    FULLCURVE_TES,
    fullcurve_signal[:, 0],
    alpha=0.5,
    color="black",
)

ax.set_ylabel("Recorded BOLD signal", fontsize=24)
ax.set_xlabel("Echo Time (ms)", fontsize=24)
ax.set_ylim(0, np.ceil(np.max(fullcurve_signal) / 1000) * 1000)
ax.set_xlim(0, np.max(FULLCURVE_TES))
ax.tick_params(axis="both", which="major", labelsize=14)
fig.tight_layout()
glue("fig_signal_decay_single-echo", fig, display=False)
../_images/Signal_Decay_10_1.png
../_images/Signal_Decay_10_0.png

Fig. 2 BOLD signal decay for a standard single-echo scan

Plot BOLD signal decay and BOLD contrast

fullcurve_signal_active = predict_bold_signal(FULLCURVE_TES, [MEAN_S0], [40])
fullcurve_signal_inactive = predict_bold_signal(FULLCURVE_TES, [MEAN_S0], [20])

fig, ax = plt.subplots(figsize=(14, 7.5))
ax.plot(
    FULLCURVE_TES,
    fullcurve_signal_active[:, 0],
    alpha=0.5,
    color="red",
    label="High Activity",
)
ax.plot(
    FULLCURVE_TES,
    fullcurve_signal_inactive[:, 0],
    alpha=0.5,
    color="blue",
    label="Low Activity",
)

ax.set_ylabel("Recorded BOLD signal", fontsize=24)
ax.set_xlabel("Echo Time (ms)", fontsize=24)
ax.set_ylim(0, np.ceil(np.max(fullcurve_signal) / 1000) * 1000)
ax.set_xlim(0, np.max(FULLCURVE_TES))
ax.legend(fontsize=20)
ax.tick_params(axis="both", which="major", labelsize=14)
fig.tight_layout()
glue("fig_signal_decay_contrast", fig, display=False)
../_images/Signal_Decay_12_1.png
../_images/Signal_Decay_12_0.png

Fig. 3 BOLD signal decay and BOLD contrast

Plot single-echo data resulting from \(S_{0}\) and \(T_{2}^{*}\) fluctuations

This shows how fMRI data fluctuates over time.

Plot single-echo data and the curve resulting from \(S_{0}\) and \(T_{2}^{*}\) fluctuations

This shows how single-echo data is a sample from a signal decay curve.

fullcurve_signal = predict_bold_signal(FULLCURVE_TES, s0_ts, t2s_ts)
singleecho_signal = fullcurve_signal[SINGLEECHO_TE, :]

# gif code based on https://www.geeksforgeeks.org/create-an-animated-gif-using-python-matplotlib/
fig, axes = plt.subplots(nrows=2, figsize=(14, 10), gridspec_kw={"height_ratios": [1, 3]})

# Set constant params for figure
axes[0].set_ylabel("Signal", fontsize=24)
axes[0].set_xlabel("Volume", fontsize=24)
axes[0].set_xlim(0, N_VOLS - 1)
axes[0].tick_params(axis="both", which="major", labelsize=14)
axes[1].set_ylabel("Signal", fontsize=24)
axes[1].set_xlabel("Echo Time (ms)", fontsize=24)
axes[1].set_xticks(MULTIECHO_TES)
axes[1].set_ylim(0, np.ceil(np.max(fullcurve_signal) / 1000) * 1000)
axes[1].set_xlim(0, np.max(FULLCURVE_TES))
axes[1].tick_params(axis="both", which="major", labelsize=14)
fig.tight_layout()

ax0_line_plot = axes[0].plot(
    singleecho_signal[0, :],
    color="black",
    zorder=0
)
ax0_scatter_plot = axes[0].scatter(
    0,
    singleecho_signal[:, 0],
    color="orange",
    s=150,
    label="Single-Echo Signal",
    zorder=1,
)
ax1_line_plot = axes[1].plot(
    FULLCURVE_TES,
    fullcurve_signal[:, 0],
    alpha=0.5,
    color="black",
    zorder=0,
)[0]
ax1_scatter_plot = axes[1].scatter(
    SINGLEECHO_TE,
    singleecho_signal[:, 0],
    color="orange",
    s=150,
    alpha=1.,
    label="Single-Echo Signal",
    zorder=1,
)

def AnimationFunction(frame):
    """Function takes frame as an input."""
    ax0_scatter_plot.set_offsets(
        np.column_stack((
            frame,
            singleecho_signal[:, frame],
        ))
    )

    ax1_line_plot.set_data((
        FULLCURVE_TES,
        fullcurve_signal[:, frame],
    ))
    ax1_scatter_plot.set_offsets(
        np.column_stack((
            SINGLEECHO_TE,
            singleecho_signal[:, frame],
        ))
    )

anim_created = FuncAnimation(fig, AnimationFunction, frames=N_VOLS, interval=100)
html = display.HTML(anim_created.to_jshtml())
glue("fig_signal_decay", html, display=False)
../_images/Signal_Decay_14_1.png

Fig. 4 Single-echo data and the curve resulting from \(S_{0}\) and \(T_{2}^{*}\) fluctuations

Plot single-echo data, the curve, and the \(S_{0}\) and \(T_{2}^{*}\) values resulting from \(S_{0}\) and \(T_{2}^{*}\) fluctuations

This shows how changes in fMRI data can be driven by both S0 and T2* fluctuations.

fullcurve_signal = predict_bold_signal(FULLCURVE_TES, s0_ts, t2s_ts)
singleecho_signal = fullcurve_signal[SINGLEECHO_TE, :]

fig, axes = plt.subplots(nrows=2, figsize=(14, 10), gridspec_kw={"height_ratios": [1, 3]})

t2s_value = predict_bold_signal(
    np.array([t2s_ts[0]]),
    np.array([s0_ts[0]]),
    np.array([t2s_ts[0]])
)[0]

ax0_line_plot = axes[0].plot(
    singleecho_signal[0, :],
    color="black",
    zorder=0
)
ax0_scatter_plot = axes[0].scatter(
    0,
    singleecho_signal[:, 0],
    color="orange",
    s=150,
    label="Single-Echo Signal",
    zorder=1,
)
ax1_line_plot = axes[1].plot(
    FULLCURVE_TES,
    fullcurve_signal[:, 0],
    alpha=0.5,
    color="black",
    zorder=0,
)[0]
ax1_scatter_plot = axes[1].scatter(
    SINGLEECHO_TE,
    singleecho_signal[:, 0],
    color="orange",
    s=150,
    alpha=1.,
    label="Single-Echo Signal",
    zorder=1,
)
ax1_t2s_scatter_plot = axes[1].scatter(
    t2s_ts[0],
    t2s_value,
    marker="*",
    color="blue",
    s=500,
    alpha=0.5,
    label="$T_{2}^{*}$",
    zorder=1,
)
ax1_s0_scatter_plot = axes[1].scatter(
    0,
    s0_ts[0],
    marker="*",
    color="red",
    s=500,
    alpha=0.5,
    label="$S_{0}$",
    clip_on=False,
    zorder=2,
)

# Set constant params for figure
axes[0].set_ylabel("Signal", fontsize=24)
axes[0].set_xlabel("Volume", fontsize=24)
axes[0].set_xlim(0, N_VOLS - 1)
axes[0].tick_params(axis="both", which="major", labelsize=14)
axes[1].legend(loc="upper right", fontsize=20)
axes[1].set_ylabel("Signal", fontsize=24)
axes[1].set_xlabel("Echo Time (ms)", fontsize=24)
axes[1].set_xticks(MULTIECHO_TES)
axes[1].set_ylim(0, np.ceil(np.max(fullcurve_signal) / 1000) * 1000)
axes[1].set_xlim(0, np.max(FULLCURVE_TES))
axes[1].tick_params(axis="both", which="major", labelsize=14)
fig.tight_layout()


def AnimationFunction(frame):
    """Function takes frame as an input."""
    t2s_value = predict_bold_signal(
        np.array([t2s_ts[frame]]),
        np.array([s0_ts[frame]]),
        np.array([t2s_ts[frame]])
    )[0]

    ax0_scatter_plot.set_offsets(
        np.column_stack((
            frame,
            singleecho_signal[:, frame],
        ))
    )

    ax1_line_plot.set_data((
        FULLCURVE_TES,
        fullcurve_signal[:, frame],
    ))
    ax1_scatter_plot.set_offsets(
        np.column_stack((
            SINGLEECHO_TE,
            singleecho_signal[:, frame],
        ))
    )
    ax1_t2s_scatter_plot.set_offsets(
        np.column_stack((
            t2s_ts[frame],
            t2s_value,
        ))
    )
    ax1_s0_scatter_plot.set_offsets(
        np.column_stack((
            0,
            s0_ts[frame],
        ))
    )

anim_created = FuncAnimation(fig, AnimationFunction, frames=N_VOLS, interval=100)
html = display.HTML(anim_created.to_jshtml())
glue("fig_signal_decay2", html, display=False)
../_images/Signal_Decay_16_1.png

Fig. 5 Single-echo data, the curve, and the \(S_{0}\) and \(T_{2}^{*}\) values resulting from \(S_{0}\) and \(T_{2}^{*}\) fluctuations

Plot \(S_{0}\) and \(T_{2}^{*}\) fluctuations

This shows how fluctuations in S0 and T2* produce different patterns in the full signal decay curves.

s0based_fullcurve_signal = predict_bold_signal(FULLCURVE_TES, s0_ts, np.full(N_VOLS, MEAN_T2S))
t2sbased_fullcurve_signal = predict_bold_signal(FULLCURVE_TES, np.full(N_VOLS, MEAN_S0), t2s_ts)

fig, axes = plt.subplots(nrows=2, figsize=(14, 10), gridspec_kw={"height_ratios": [1, 3]})

t2s_value = predict_bold_signal(
    np.array([t2s_ts[0]]),
    np.array([s0_ts[0]]),
    np.array([t2s_ts[0]])
)[0]

ax0_line_plot = axes[0].plot(
    ts,
    color="black",
    zorder=0
)
ax0_scatter_plot = axes[0].scatter(
    0,
    ts[0],
    color="purple",
    s=150,
    label="Single-Echo Signal",
    zorder=1,
)
ax1_s0_line_plot = axes[1].plot(
    FULLCURVE_TES,
    s0based_fullcurve_signal[:, 0],
    alpha=0.5,
    color="red",
    label="$S_0$-Driven Decay Curve",
)[0]
ax1_t2s_line_plot = axes[1].plot(
    FULLCURVE_TES,
    t2sbased_fullcurve_signal[:, 0],
    alpha=0.5,
    color="blue",
    label="$T_{2}^{*}$-Driven Decay Curve",
)[0]

# Set constant params for figure
axes[0].set_ylabel("$S_0$/$T_{2}^{*}$", fontsize=24)
axes[0].set_xlabel("Volume", fontsize=24)
axes[0].set_xlim(0, N_VOLS - 1)
axes[0].tick_params(axis="both", which="major", labelsize=14)
axes[1].legend(loc="upper right", fontsize=20)
axes[1].set_ylabel("Signal", fontsize=24)
axes[1].set_xlabel("Echo Time (ms)", fontsize=24)
axes[1].set_xticks(MULTIECHO_TES)
axes[1].set_ylim(0, np.ceil(np.max(s0based_fullcurve_signal) / 1000) * 1000)
axes[1].set_xlim(0, np.max(FULLCURVE_TES))
axes[1].tick_params(axis="both", which="major", labelsize=14)
fig.tight_layout()


def AnimationFunction(frame):
    """Function takes frame as an input."""
    ax0_scatter_plot.set_offsets(
        np.column_stack((
            frame,
            ts[frame],
        ))
    )

    ax1_t2s_line_plot.set_data((
        FULLCURVE_TES,
        t2sbased_fullcurve_signal[:, frame],
    ))
    ax1_s0_line_plot.set_data((
        FULLCURVE_TES,
        s0based_fullcurve_signal[:, frame],
    ))

anim_created = FuncAnimation(fig, AnimationFunction, frames=N_VOLS, interval=100)
html = display.HTML(anim_created.to_jshtml())
glue("fig_signal_decay3", html, display=False)
../_images/Signal_Decay_18_1.png

Fig. 6 \(S_{0}\) and \(T_{2}^{*}\) fluctuations

Plot \(S_{0}\) and \(T_{2}^{*}\) fluctuations and resulting single-echo data

This shows how single-echo data, on its own, cannot distinguish between S0 and T2* fluctuations.

s0based_fullcurve_signal = predict_bold_signal(FULLCURVE_TES, s0_ts, np.full(N_VOLS, MEAN_T2S))
t2sbased_fullcurve_signal = predict_bold_signal(FULLCURVE_TES, np.full(N_VOLS, MEAN_S0), t2s_ts)
s0based_singleecho_signal = s0based_fullcurve_signal[SINGLEECHO_TE, :]
t2sbased_singleecho_signal = t2sbased_fullcurve_signal[SINGLEECHO_TE, :]

fig, axes = plt.subplots(nrows=2, figsize=(14, 10), gridspec_kw={"height_ratios": [1, 3]})

ax0_line_plot = axes[0].plot(
    ts,
    color="black",
    zorder=0
)
ax0_scatter_plot = axes[0].scatter(
    0,
    ts[0],
    color="purple",
    s=150,
    zorder=1,
)
ax1_s0_line_plot = axes[1].plot(
    FULLCURVE_TES,
    s0based_fullcurve_signal[:, 0],
    color="red",
    alpha=0.5,
    label="$S_0$-Driven Decay Curve",
    zorder=0,
)[0]
ax1_t2s_line_plot = axes[1].plot(
    FULLCURVE_TES,
    t2sbased_fullcurve_signal[:, 0],
    color="blue",
    alpha=0.5,
    label="$T_{2}^{*}$-Driven Decay Curve",
    zorder=1,
)[0]
ax1_s0_scatter_plot = axes[1].scatter(
    SINGLEECHO_TE,
    s0based_singleecho_signal[:, 0],
    color="red",
    s=150,
    alpha=1.,
    label="$S_0$-Driven Single-Echo Signal",
    zorder=2,
)
ax1_t2s_scatter_plot = axes[1].scatter(
    SINGLEECHO_TE,
    t2sbased_singleecho_signal[:, 0],
    color="blue",
    s=150,
    alpha=1.,
    label="$T_{2}^{*}$-Driven Single-Echo Signal",
    zorder=3,
)

# Set constant params for figure
axes[0].set_ylabel("$S_0$/$T_{2}^{*}$", fontsize=24)
axes[0].set_xlabel("Volume", fontsize=24)
axes[0].set_xlim(0, N_VOLS - 1)
axes[0].tick_params(axis="both", which="major", labelsize=14)
axes[1].legend(loc="upper right", fontsize=20)
axes[1].set_ylabel("Signal", fontsize=24)
axes[1].set_xlabel("Echo Time (ms)", fontsize=24)
axes[1].set_xticks(MULTIECHO_TES)
axes[1].set_ylim(0, np.ceil(np.max(s0based_fullcurve_signal) / 1000) * 1000)
axes[1].set_xlim(0, np.max(FULLCURVE_TES))
axes[1].tick_params(axis="both", which="major", labelsize=14)
fig.tight_layout()


def AnimationFunction(frame):
    """Function takes frame as an input."""
    ax0_scatter_plot.set_offsets(
        np.column_stack((
            frame,
            ts[frame],
        ))
    )

    ax1_s0_line_plot.set_data((
        FULLCURVE_TES,
        s0based_fullcurve_signal[:, frame],
    ))
    ax1_t2s_line_plot.set_data((
        FULLCURVE_TES,
        t2sbased_fullcurve_signal[:, frame],
    ))
    ax1_s0_scatter_plot.set_offsets(
        np.column_stack((
            SINGLEECHO_TE,
            s0based_singleecho_signal[:, frame],
        ))
    )
    ax1_t2s_scatter_plot.set_offsets(
        np.column_stack((
            SINGLEECHO_TE,
            t2sbased_singleecho_signal[:, frame],
        ))
    )

anim_created = FuncAnimation(fig, AnimationFunction, frames=N_VOLS, interval=100)
html = display.HTML(anim_created.to_jshtml())
glue("fig_signal_decay4", html, display=False)
../_images/Signal_Decay_20_1.png

Fig. 7 \(S_{0}\) and \(T_{2}^{*}\) fluctuations and resulting single-echo data

Plot \(S_{0}\) and \(T_{2}^{*}\) fluctuations and resulting multi-echo data

This shows how S0 and T2* fluctuations produce different patterns in multi-echo data.

s0based_fullcurve_signal = predict_bold_signal(FULLCURVE_TES, s0_ts, np.full(N_VOLS, MEAN_T2S))
t2sbased_fullcurve_signal = predict_bold_signal(FULLCURVE_TES, np.full(N_VOLS, MEAN_S0), t2s_ts)
s0based_multiecho_signal = s0based_fullcurve_signal[MULTIECHO_TES, :]
t2sbased_multiecho_signal = t2sbased_fullcurve_signal[MULTIECHO_TES, :]

fig, axes = plt.subplots(nrows=2, figsize=(14, 10), gridspec_kw={"height_ratios": [1, 3]})

ax0_line_plot = axes[0].plot(
    ts,
    color="black",
    zorder=0,
)[0]
ax0_scatter_plot = axes[0].scatter(
    0,
    ts[0],
    color="purple",
    s=150,
    zorder=1,
)
ax1_s0_line_plot = axes[1].plot(
    FULLCURVE_TES,
    s0based_fullcurve_signal[:, 0],
    color="red",
    alpha=0.5,
    label="$S_0$-Driven Decay Curve",
    zorder=0,
)[0]
ax1_t2s_line_plot = axes[1].plot(
    FULLCURVE_TES,
    t2sbased_fullcurve_signal[:, 0],
    color="blue",
    alpha=0.5,
    label="$T_{2}^{*}$-Driven Decay Curve",
    zorder=1,
)[0]
ax1_s0_scatter_plot = axes[1].scatter(
    MULTIECHO_TES,
    s0based_multiecho_signal[:, 0],
    color="red",
    s=150,
    alpha=1.,
    label="$S_0$-Driven Multi-Echo Signal",
    zorder=2,
)
ax1_t2s_scatter_plot = axes[1].scatter(
    MULTIECHO_TES,
    t2sbased_multiecho_signal[:, 0],
    color="blue",
    s=150,
    alpha=1.,
    label="$T_{2}^{*}$-Driven Multi-Echo Signal",
    zorder=3,
)

axes[0].set_ylabel("$S_0$/$T_{2}^{*}$", fontsize=24)
axes[0].set_xlabel("Volume", fontsize=24)
axes[0].set_xlim(0, N_VOLS - 1)
axes[0].tick_params(axis="both", which="major", labelsize=14)
axes[1].legend(loc="upper right", fontsize=20)
axes[1].set_ylabel("Signal", fontsize=24)
axes[1].set_xlabel("Echo Time (ms)", fontsize=24)
axes[1].set_xticks(MULTIECHO_TES)
axes[1].set_ylim(0, np.ceil(np.max(s0based_fullcurve_signal) / 1000) * 1000)
axes[1].set_xlim(0, np.max(FULLCURVE_TES))
axes[1].tick_params(axis="both", which="major", labelsize=14)
fig.tight_layout()


def AnimationFunction(frame):
    """Function takes frame as an input."""
    ax0_scatter_plot.set_offsets(
        np.column_stack((
            frame,
            ts[frame],
        ))
    )

    ax1_s0_line_plot.set_data((
        FULLCURVE_TES,
        s0based_fullcurve_signal[:, frame],
    ))
    ax1_t2s_line_plot.set_data((
        FULLCURVE_TES,
        t2sbased_fullcurve_signal[:, frame],
    ))
    ax1_s0_scatter_plot.set_offsets(
        np.column_stack((
            MULTIECHO_TES,
            s0based_multiecho_signal[:, frame],
        ))
    )
    ax1_t2s_scatter_plot.set_offsets(
        np.column_stack((
            MULTIECHO_TES,
            t2sbased_multiecho_signal[:, frame],
        ))
    )

anim_created = FuncAnimation(fig, AnimationFunction, frames=N_VOLS, interval=100)
html = display.HTML(anim_created.to_jshtml())
glue("fig_signal_decay5", html, display=False)
../_images/Signal_Decay_22_1.png

Fig. 8 \(S_{0}\) and \(T_{2}^{*}\) fluctuations and resulting multi-echo data

Plot \(T_{2}^{*}\) against BOLD signal from single-echo data (TE=30ms)

When the BOLD signal is driven entirely by T2* fluctuations, the BOLD signal is very similar to a scaled version of those T2* fluctuations, but there are small differences.

fig, ax = plt.subplots(figsize=(16, 4))

scalar = np.linalg.lstsq(
    t2s_ts[:, None],
    t2sbased_fullcurve_signal[SINGLEECHO_TE[0], :],
    rcond=None
)[0]
ax.plot(
    t2sbased_fullcurve_signal[SINGLEECHO_TE[0], :],
    color="orange",
    label=f"BOLD signal (TE={SINGLEECHO_TE[0]}ms)",
)
ax.plot(
    t2s_ts * scalar,
    label="$T_{2}^{*}$ (scaled to signal)",
    linestyle="--",
)
ax.set_xlim(0, N_VOLS - 1)
leg = ax.legend()
glue("fig_t2s_bold_single-echo", fig, display=False)
../_images/Signal_Decay_24_1.png
../_images/Signal_Decay_24_0.png

Fig. 9 \(T_{2}^{*}\) against BOLD signal from single-echo data (TE=30ms)

Plot \(S_{0}\) against BOLD signal from single-echo data (TE=30ms)

When the BOLD signal is driven entirely by S0 fluctuations, the BOLD signal ends up being a scaled version of the S0 fluctuations.

fig, ax = plt.subplots(figsize=(16, 4))

scalar = np.linalg.lstsq(
    s0_ts[:, None],
    s0based_fullcurve_signal[SINGLEECHO_TE[0], :],
    rcond=None
)[0]
ax.plot(
    s0based_fullcurve_signal[SINGLEECHO_TE[0], :],
    color="orange",
    label=f"BOLD signal (TE={SINGLEECHO_TE[0]}ms)",
)
ax.plot(
    s0_ts * scalar,
    label="S0 (scaled to signal)",
    linestyle="--",
)
ax.set_xlim(0, N_VOLS - 1)
leg = ax.legend()
glue("fig_s0_bold_single-echo", fig, display=False)
../_images/Signal_Decay_26_1.png
../_images/Signal_Decay_26_0.png

Fig. 10 \(S_{0}\) against BOLD signal from single-echo data (TE=30ms)